%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Empirical Application(Timber Dataset)               %
% Test: Exogenous Bidders Participation               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all

load('data.mat'); % load timber dataset
L2=size(wb2,1); % number of auctions with N=2
L3=size(wb3,1); % number of auctions with N=3
L=L2+L3;

% Covariates
x2=[ones(L2,1),apv2,vol2];
x3=[ones(L3,1),apv3,vol3];
x = [x2;x3];
p=size(x,2);

% Winning bids
wb=[wb2;wb3];

N=[2,3];
alpha=0.12:0.02:0.80; % quantile levels alpha
m=size(alpha,2);
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];

B=1000; % Bootstrap replications
toler=0.1; % Parameter of tolerance for precision of the estimation
gammai=ones(p,1); % initial guess for gamma

% Computation of the Exogenous Bidder Participation Test Statistic
for k=1:m
    % Estimation under H0
    [obj0,res0,~]=MM(psi(:,k), p, wb, x, toler, gammai);
    % residuals of the estimation under H0
    res0k(:,k)=res0; 
    % Estimation under H1
    [obj2,res2,~]=MM(psi2(:,k),p,wb2,x2,toler,gammai);
    [obj3,res3,~]=MM(psi3(:,k),p,wb3,x3,toler,gammai);
    % residuals of the estimation under H1
    res2k(:,k)=res2;
    res3k(:,k)=res3;
    M(k,:)=2*(obj0-(obj2+obj3));
end
   
   Mind=sum(M); % test statistic
   
% Random Weighting Bootstrap of the Exogenous Bidder Participation Test Statistic  
   
for b=1:B
    
    w=[mnrnd(L2,ones(L2,1)/L2)';mnrnd(L3,ones(L3,1)/L3)']; % weights
                
    for k=1:m
        % residuals under H0 weighted by w
        res0w=res0k(:,k).*w; 
        % residuals under H1 weighted by w
        res2w=res2k(:,k).*w(1:L2,:);
        res3w=res3k(:,k).*w(L2+1:L,:); 
        % objective functions under H0 weighted by w
        obj0w=sum(psi(:,k).*res0w)-sum(res0w(res0w<0)); 
        % objective functions under H1 weighted by w
        obj2w=sum(psi2(:,k).*res2w)-sum(res2w(res2w<0));
        obj3w=sum(psi3(:,k).*res3w)-sum(res3w(res3w<0));
        Mw(k,:)=2*(obj0w-(obj2w+obj3w));
        % objective functions under the random weighting method
        [obj0wb,~,~]=MMw(psi(:,k),p,wb,x,w,toler,gammai); % under H0
        [obj2wb,~,~]=MMw(psi2(:,k),p,wb2,x2,w(1:L2,:),toler,gammai); % under H1
        [obj3wb,~,~]=MMw(psi3(:,k),p,wb3,x3,w(L2+1:L,:),toler,gammai); % under H1
        Ms(k,:)=2*(obj0wb-(obj2wb+obj3wb));
    end
        
    Mstar(:,b)=sum(Ms-Mw);
    
    b    
end
    
  % Critical Values
  
    c99=quantile(Mstar',0.99);
    c95=quantile(Mstar',0.95);
    c90=quantile(Mstar',0.90);
    
  % P-value computation
  
     R=zeros(1,B);
for b=1:B
    if Mstar(:,b)>=Mind
            R(:,b)=1;
    else
            R(:,b)=0;
    end
end

    pvalue=sum(R')/B
    